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INTRODUCTION 


The expanding uses of composite materials, especially of the 
filamentary type, require the generation of a variety of data on 
material behavior so that adequate design can be achieved. Many 
test techniques have been developed for such purposes. Among these 
are the nondestructive evaluation (NDE) techniques of ultrasonics 
and acoustic emission. 

The application and interpretation of results obtained from various 
NDE techniques are strongly related to the stress wave propagation 
characteristics of the material. The theoretical prediction of the 
stress wave fields under the action of well defined excitation can 
give information regarding the proper use of the NDE techniques, and 
may provide improved precision of a material's diagnosis by 
establishing more quantitative standards. 

Via an asymptotic approximation, the stress wave field generated 
by a point source in an infinite anisotropic medium can be described 
in analogy with the problem of magnetohydrodynamic waves discussed 
by Lighthill [1] . For the present work, the particular case of 
transversely isotropic media is presented by following the same 
procedure used by Buchwald [2], where here the solution is extended 
such that the displacements are obtained and expressed in cartesian 
form. 

Here, the solution is given for the case of a glass fiber 
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reinforced epoxy unidirectional composite. The slowness and wave 
surfaces are determined as well as the displacements for points 
throughout the medium. Polar diagrams for displacement amplitudes 
are constructed illustrating the patterns of the displacement field. 

Other aspects of the problem of stress waves in homogeneous 
anisotropic media can be found in the studies of Synge [3] , Carrier 
[4] and Musgrave [5]. Synge analyzed the behavior of stress waves 
in a single layer of anisotropic material subjected to uniform stress 
at one of its surfaces and extended the study to the limiting case 
where the layer is transformed into a half space. Carrier proposed 
a solution for waves in an infinite medium using potential functions 
and integral techniques. Musgrave established the conditions to be 
satisfied for the existence of the general equation of motion 
when a plane wave solution is adopted and described the geometry of 
the spreading disturbances by means of the wave surfaces and 
slowness surfaces. 

The work on layered media was carefully explored but was not 
pursued for this research. The primary reason for this is that 
general anisotropy is a more convenient approach for representing 
a composite because the effective modulus theory can be used to 
model the material as homogeneous [6], It is assumed that the 
material behaves as homogeneous since the layer thicknesses (fiber 
diameter for filamentary composites) are small as compared with the 
wavelengths. If this requirement is satisfied, the effective 
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modulus theory appears to describe the behavior of stress waves in 
laminated composites better than the laminated media theory [7,8], 
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EQUATION OF MOTION AND ITS SOLUTION 


A schematic of the system under investigation is shown in Fig. 1. 
An infinite transversely isotropic medium is subjected to a point 
force source which generates stress waves that propagate throughout 
the medium. A cartesian coordinate system x, y and z is adopted such 
that the x-y plane is the isotropic plane in the medium. In a 
unidirectional composite, z would be the fiber direction. The point 
source is located at the origin and oscillates in the x direction as 
shown. The displacements at any point P in the medium due to a 
steady-state sinusoidal point force excitation are analyzed. 

Force-Dynamic Equations 

The force-dynamic equations for an infinite medium can be derived 
from the equilibrium conditions of a differential volumetric element 
and are given as [9] 


Pu. = T . + p F. 
J ,tt ij ,i j 


( 1 ) 


where the implied summation is over the repeated subscripts, p is the 

mass density, the indices i and j have values x, y and z in a 

rectangular cartesian coordinate system, uj is the displacement in the 

i direction, x . . is the stress tensor and F. is the body force (force 
ij J 

per unit mass) in the j direction. The indices following the comma 
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refer to partial derivatives such that u^ tt represents the second 
derivative of the jth component of displacement with respect to t, 
the time variable. 

For the point sinusoidal force shown in Fig. 1, the body 
forces can be written in exponential form as 

F = F 6 (x) 6 (y) 6 (z) e ia)t and F = F =■ 0 

x o y z 

( 2 ) 

where i = (-1), oj is the radian frequency of the harmonic excitation, 

6( ) is the Dirac delta function and F is the force amplitude. The 

o 

actual oscillatory force is taken as the real part of the expression 
for F^ in eqn. (2). 

Constitutive Relations 

The constitutive relations based on Hooke’s law are [10] 

x . . = C. . e (3) 

upq pq 

where C . are the elastic coefficients that constitute the stiffness 

ijpq 

matrix and e is the strain tensor that can expressed in terms of 

Pq 

displacements using the engineering definitions for strains as 


epq 



p.q 


+ u 

q,p 


) 


(4) 


where the suffixes each represent the x, y and z directions. The 
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elastic coefficients satisfy the symmetry conditions such that [11] 


c ijpq C jipq c ijqp = c pqij 


(5) 


The four indices notation for the elastic coefficients can be 
modified to the two indices compact form [9]. In this case, the 
stress tensor is written as a vector 

T 1 “ t T l’ t 2 , t 3 , t 4 , t 5 , Tg] (6) 


where the equivalence in the cartesian system is = T xx » T 2 3 T yy’ 

t 3 = T zz’ t 4 “ T X z’ T 5 = T yz’ T 6 = T xy* The same corres P onden ce 
is applicable for the strains. 

For a transversely isotropic medium as shown in Fig. 1, where 
the x - y plane is the isotropic plane, z is the symmetry axis 
of the medium. It can be shown that for transversely isotropic 
materials, there are five independent elastic coefficients [11]. 

The nonzero elastic coefficients according to the definition of 
the stress components in eqn. (6) are 


C 11 “ C 22’ C 33* C 12* C 13 


C 23 * C 44 


C 55 and C g6 - 2 (C u - C 12 > 


(7) 


Constraints 


The constraints are introduced by the application of the 
"radiation condition", Since the medium is infinite, the constraint 
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is that no waves originate at infinity or are reflected at infinity; 
that is, the only source of disturbance is due to the point force 
located at the origin of the coordinate system. From this 
assumption, the amplitudes of the displacements must decrease with 
distance from the source. In the limit, as the distance approaches 
infinity, the displacement amplitudes must vanish. 


Equations of Motion 

The equations of motion can be represented in terms of the 
displacements by substituting eqns. (4) and (5) into eqn. (3) and 
then substituting the result into eqn. (1) to give [10] 


pu 


= C 


u 


j,tt ijpq p,qi 


+ p 


F. 

J 


( 8 ) 


where the symmetry condition given in eqn. (5) is utilized. The 
forcing function is as given in eqn. (2) . For a transversely 
isotropic medium, the elastic coefficients obey eqn. (7) . 

Solution of Equations of Motion 

The equations of motion are given in eqn. (8) and must be 
solved in accordance with the imposed constraints and the applied 
forces described earlier. 

Eqn. (8) can be rewritten using u, v and w for the 
displacements in the x, y and z directions, respectively, instead 
of u x> u^ and u^; and using x, y and z for the subscripts denoting 
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derivatives. This gives three equations in u, v and w. Then 

these equations can be differentiated with respect to the space 

and time variables to give terms in u, , u, . v, , v. 

a ’ttx’ ’tty’ ttx ’tty 

and w, which can be combined by addition and subtraction to 

UuZ 

give [2] 


A, 

r, 

A, 


tt 

tt 

tt 


a.V. A + a_ A, + F 
4 1 5 zz x,y 


= a 3 A ’zz + a 5 A 1 F + a ’ r ’ 


2 ’zz 


- a.V.r + a_ A, + a. V.A + F 
3 1 5 ’zz 1 1 x,x 


(9) 

( 10 ) 

( 11 ) 


where the new variables are 


A = 





r 


= W, 


A = 


x 


+ V, 


( 12 ) 


A represents the z component of rotation, T is the strain in the 
z direction and (T + A) is the dilatation, all in accordance with 
standard definitions in elasticity. Physically, the fact that 
a separate equation for the z rotation can be written, namely eqn. 
(9) , indicates that the propagation of these rotational waves is 
independent of the other possible modes. The same is not true 
for the displacement functions T and A, since eqns. (10) and (11) 
cannot be decoupled. The a's are constants given by 
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S 3 G 44 4€ 13 C 11~ C 12 C 44 

a l = p * 2 " p ’ 3 ~ p * 4 “ 2 P * a 5 “ P 

(13) 


and is the laplacian operator in two dimensions such that 


V 


2 

1 



(14) 


If eqn. (2) is substituted into eqns. (9), (10) and (11), the 

asymptotic solution can be obtained by using Fourier integrals and 

the method of stationary phase [12]. Basically, eqns. (9), (10) 

and (11) are written in terms of Fourier integrals. In doing this, 

plane wave solutions are assumed with associated wave number vectors 

that are represented by their components a along the x direction, 

B along the y direction and y along the z direction. Then, the 

values of the Fourier integrals are approximated by their residues 

such that the stationary points are determined by the singularities 

of the integrands. In this approximation, terms involving the 
2 

factor 1/r (where r is the radial distance from the origin to P the 
point of interest) are neglected. Thus, the expressions obtained 
are adequate for the calculation of the displacements of points 
in the thus defined far field. This leads to the following results [2] 

^ i S F 

A = — r-py exp [i(ax + By + Yz - “t)] (15) 

2rr|VG| (|K G |) i/Z 
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where the symbol | | denote "the magnitude of". Now a, 6 and y are 

the wave number components defining the singularities of the integrands 
in the integration process for which the residues are calculated. 

Eqn. (15) represents the residue corresponding to a single singularity 
point. For the eqns. (16) and (17) the wave number components are 
followed by the index n indicating that it is possible to have more 
than one singularity point and in this case the expressions for A and 
F represent the sum of the residues corresponding to each of the 
singularity points. N is the total number of singularity points over 
which the summation must be performed. Also 


2 2 2 2 
a^,(a + 8 ; + a^y - cj 


|VG| = [Gf a + Gfg + G, y ] (19) 

2 222 22 222222 
H = [a 5 y>a 1 (a Z +6 Z )-a) Z ][a 5 (a^+B )+a 2 y Z -(o Z ]-a3Y^(a Z +B ) 

( 20 ) 
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( 21 ) 


|7H| = [H^ + H? g + H? y ] 


1/2 


K_ = 


G, 00 [G? C, -2G, G, C, + G 2 , G, ] 
SB a YY a y ’ay ’y ’aa 

< + G V 2 


( 22 ) 




H, d An 2 , H, -2H, H, H, +H^ H, ] 
BB a YY a Y ay ’y aa 

2 2 2 
(H» a + Hf y ) Z 


(23) 


2 2 . 2 , 2 
r = x + y + z 


(24) 


The values and are commonly known as gaussian curvatures [1] , 
and the combination of |vh| and ( | K^|) after simplification [2] 
is usually expressed by 


A = 


n ItoKIkJ) 


1/2 


2 2 
h; + h; 

’a y 


1 1/2 


H, q 0 (H^ H, -2H, H, H, +H? H, ) 
’BB a yy a y ay y aa 


(25) 


A is a phase coefficient that is determined by the geometric properties 
n 

of the H = 0 surface as follows: 

a) A = + 1 if K > 0 and VH is in the +r direction; 

n H — 

b) A = - 1 if K > 0 and VH is in the -r direction; 

n h — 
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c) A = + 1 if K < 0 and H = 0 is convex in the VH direction; 

n n — 

and 

d) A = - 1 if K < 0 and H = 0 is concave in the VH direction. 

n H — 

The values for A^ can be determined only by numerical calculations. 

As mentioned before, eqns. (15), (16) and (17) are residues or 
sums of residues resulting from the integration process. In wave 
number space, G = 0 and H = 0 are surfaces containing all possible 
values of wave numbers that are able to generate such residues, 
considering the elastic properties of the material. Each point 
(a, 3, y) on the surfaces can then be associated with a particular 
wave train in a particular direction which effectively contributes 
to the displacement functions. It is important to observe that the 
H = 0 surface is in reality a quartic of two sheets [2] ; so H = 0 
defines two surfaces. 

According to the method of stationary phase, the displacements 

at a point P in the medium are calculated by summing the contributions 

of plane waves passing through such a point [1] . These waves and the 

corresponding wave numbers can be found by locating the point 

(a, 3, y) or points (a , 3 , Y ) on the surfaces G = 0 and H = 0, 

n n n 

respectively, where the normal is parallel to the direction OP (in 
Fig. 1), which corresponds to finding the singularity points that 
generate the residues for the specific point P. 


- 12 - 



Note that normals to the G = 0 or H =* 0 surfaces at more than 


one point may lie in the same direction (parallel to OP) . In such 
a direction waves of different settings (direction of wave number 
vector) and spacings (magnitude of wave number vector) can be super- 
posed and for this reason the summation signs appear in eqns. (16) 
and (17) . It must be noted that this remark does not apply if the 
points are simply symmetric such as (a, 6, y) and (-a, -8, -y) which 
correspond to waves of the same setting and spacing. For the solution 
satisfying the radiation condition, only one out of each such pair 
of points is selected. Geometrically, the surface G =» 0 for a 
transversely isotropic medium is an ellipsoid and as a consequence 
there is only one point on the surface where the normal is parallel 
to OP. 

With the use of eqns. (18) through (25), the expressions given 
in eqns. (15), (16) and (17) for the displacement functions can be 
rewritten as 


2 2 2 2 2 
* 6 F Q [a^ (a + g Z ) + a* Y ] 

= " ' 77172 

4irr w (a 5 Q) a^ 


exp [i(ax + 6y + yz - wt)] 


(26) 


. i N 

A = T — y A i X F a [a,.(a^+8^) + a_y 4 - u> 4 ] 
2jrr n n o n 5 n n' 2 n 

n=l 


2 , „ 2 . 


exp [Ka^x + B ff y + y^z - ut) ] 


(27) 
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N 


l A i X F a Y exp [i (a x + £ y + y z - cot)] 


2^r " n non n 

n=l 


n n n 


(28) 


where 


2 2 2 2 2 
Q = a 4 (a + 0 ) + a* y 


(29) 


The displacements u, v and w can be determined by assuming forms 
compatible with eqns. (26), (27) and (28) such that eqns. (12) are 
satisfied. Since the displacement functions are of exponential 
form, identical exponential forms are chosen to represent the 
displacements as 


u =■ k x A + k 2 & + k 3 (30) 

v = k 4 A + k 5 A + kg (31) 

w =« k ? T + k 8 (32) 

where through kg are constants to be determined. If the 

derivatives of eqns. (30), (31) and (32) are calculated and terms 
involving the factor 1/r 2 are neglected, the substitution of the 
resulting expressions for the derivatives into eqns. (12) gives the 
values of these constants as 
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-B 

2 2 
i(g + a ) 


a 

k, = , 2. 
4 i(g + a ) 


a 

n 


*<“n + 6 n> 


( 33 ) 


k c = 


3 n 

«“ n 2 + e n 2 > 


(34) 



The constants k_ , k, and k Q are determined using the radiation 

JO o 

condition presented earlier. Since there are no sources in the medium 
except at the origin, the displacements at infinity must vanish. 

This is possible only if the constants k^, kg and kg are zero; so 


k 3 = k 6 = k 8 = 0 


(36) 


Therefore the final expressions for the displacements are 


„2 _ r 2 , 2 , „ 2 . L 22 . 

^ 3 F q [a 4 (a + 6 ) + a 5 y ] 

u = — Tjo — exp [i (ax + By + yz - tot)] 

(3 + a ) 4Trr to (a Q) 7 a 4 


N 

l 


a 


n 


2irr S ( 2 j.a 2 \ 
n=l (a +6 ) 
n n 


A. X F [a (a 2 +3 2 )+a y 2 -<o 2 ] 
n n o 5 n n 2 n 


exp[i(a x+Sy + yz - tot) ] 
r n n' n 


(37) 
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2 2 2 2 2 

* “6 V a 4 (o +B } + vi 

V = - — ~ — j 77-5 — exp [i (ax + 8y + y - wt) ] 

(a +8 ;47rr co(a 5 Q) X/ 2 a 4 2 


+ 


1 

2Trr 


N 


I 

n=l 


CtB 

n n 


2 2 
(a +8 ) 
v n n 


A X F [a (a^+6^) + a. y^ 
n n o 5 n n 2 n 


2 


] 


exp [i(ax + 8y + y z - cot)] 
r n n^ n 

(38) 


'V/ 

w = 


2:rr 


N 

l 

n=l 


n 


X F 
n o 


Y n a 3 


exp 


tI( v 


8 y + 
n-' 


Y n Z 


- (Ot)] 


(39) 


Observe that a different procedure could have been followed 
to obtain the expressions for the displacements. Recall, that the 
expressions in eqns. (12) relate the displacements u, v and w with 
the displacement functions A, T and A. Then, from the expressions 
in eqns. (26), (27) and (28), the system represented by eqn. (12) 
could have been solved by direct integration. This procedure 
is mathematically more complex than the method of assuming 
exponential forms for u, v and w and determining constants as was 
done here. 
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PHASE VELOCITY AND WAVE SURFACE 

In order to understand the propagation behavior of waves in 
anisotropic media, some concepts concerning the plane wave propagation 
characteristics are reviewed below. For this, consider the general 
equation of motion for an anisotropic medium when no body forces are 
present. From eqn. (8) withF^ ■ 0 

(40) 

(41) 


C.. u . - pu, = 0 

ijpq p.qi j.tt 

Assume a plane wave front of the form [5] 

2 

(j> =» <f>[n • x - |cj t] 


where <j> represents a characteristic with unit normal n^ and scalar 
phase velocity |cj along the normal direction. Observe that the 
characteristics represent the wave front in time and space. By the 
theory of characteristics, the displacements are functions of 4> such 
that the spatial dependence of u is as n^ « x . 

For the nontrivial solution (that is, for the propagation of 
the front) the "characteristic condition" must be applied. This 
condition can be written as [5] 


C. . $, , <f>, — D <P , <5 . = 0 

ijpq i q t ip' 


or substituting the derivatives according eqn. (41) into eqn. (42) 
gives 
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(43) 


C. . n. n - 

1 ijpq i q 


»l =! 2 « lp l 


= 0 


where n. and n are the projections (direction cosines) of n in the 
1 q — r 

directions i and q, respectively, and is the Kronecker delta. 

The solution of eqn. (43) is an eigenvalue problem for the phase 

velocity |cj for any specified direction n^. In general, there are 

three velocities (that is, eigenvalues), distinct in magnitude for each 

n . 

— r 

Now suppose that the actual point source introduces the excitation 
into the medium. If composition of plane wave fronts is used to 
represent the actual disturbance front, according to Huygen’s principle 
[5] , the front can be traced out by the envelope of plane wave fronts 
as if they were generated at the origin in all possible directions, 
at time t = 0. So the normals n^ of the component plane fronts must 
assume all possible directions in order to sweep the entire space 
around the source. For each direction n a new determinant of the 
form shown in eqn. (42) is generated and three new eigenvelocities can 
be determined. Polar diagrams can be constructed for each of the 
three velocities as functions of the corresponding normal direction. 

The surfaces expressed in polar representation (diagrams) are called 
the velocity surfaces. Moreover, the envelope of plane wave fronts 
at unit 'time is called the wave surface [3]. There is a correspondence 
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between the velocity surface and the wave surface such that if three 
velocity surfaces exist, the same number of wave surfaces can be 
generated. The existence of three surfaces indicates that there are 
three possible modes of propagation, each mode with its particular 
phase velocities. If just one mode were to exist, the physical meaning 
of the wave surface could simply be stated as the boundary between 
the disturbed and undisturbed regions in the medium at unit time after 
the source is set into action. 

Consider now the eigenvectors corresponding to the eigenvelocities. 
The symmetry of the forms in eqn. (42) indicates that the eigenvectors 
(that is, displacement modes) corresponding to each of the three 
eigenvalues are mutually orthogonal but in general none of them is 
parallel to the wave front normal direction n [2,5]. If one of the 
eigenvectors is coincident with the wave normal direction i^., the 
corresponding displacement mode is purely longitudinal and the above - 
mentioned orthogonality ensures that the remaining displacement modes 
are purely transverse. This condition is fulfilled for all wave 
propagation directions in isotropic media where each of the wave 
surfaces can be associated with a pure propagation mode, namely, 
one longitudinal mode and two transverse modes with perpendicular 
polarization planes. 

For the transversely isotropic medium as opposed to the isotropic 
medium, the eigenvectors are not coincident with the wave normals 
in all directions, so plane wave fronts travel along directions 
oblique to their normal n [5] . This means that the modes of 
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propagation are combinations of longitudinal and transverse modes. 
Exceptions occur for certain directions. These directions are: 

(1) the principal direction perpendicular to the isotropic plane and 

(2) any direction contained in the isotropic plane. The values of 
the velocities for these directions, according to the convention 
adopted for the coordinate axes, are given in Table 1. 
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SLOWNESS SURFACE 


The condition given in eqn. (43) for the existence of the solution 
of the equation of motion given in eqn. (40) can also be represented in 
terms of slownesses. The slowness _s is generally defined as a vector 
whose magnitude is the inverse of the velocity magnitude |cj , and whose 
direction is the same direction of £, or 

|sj * 1/ 1 cj (44) 


For a system with three different velocities there are three such 
slowness vectors. 

According to the form of the argument of the wave front function 
in eqn. (41) and the definition in eqn. (44) , the component of £ in 
the direction i, s^, can be written as 

S i = n i^ I— I ( 4 5) 

Substituting eqn. (45) into eqn. (43) gives 


C. . s. s. - p 6 . = 0 

ijpq i 3 IP 


(46) 


The solution for the slownesses is the same as the solution for 
the velocities. As before, if each direction in space is considered, 
three surfaces can be traced out which are the reciprocals of the 
velocity surfaces [3]. The resulting surfaces are called slowness 
surfaces. Each radius vector from the origin to the surface on the 
slowness surface has its corresponding inverse on the velocity surface. 
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The slowness surface is then an alternative way of geometrically 
representing the propagation characteristics of materials. The solution 
for slowness is often preferred to the solution for velocities since 
the former involves simpler algebraic transformations. This argument 
is supported by the fact that the determinant from eqn. (43) leads to 
an expression of twelfth order in the velocities as opposed to the 
determinant from eqn. (46) which is a sixth order equation in s^[5]. 

Recall, the case under study is that of a transversely isotropic 
medium subjected to an oscillatory point source. The motion is 
represented by eqns. (9), (10) and (11) in terms of the variables A, 

T, and A. The slowness surfaces can be identified when the differential 
equations (9), (10) and (11) are written in terms of their Fourier 
transforms [1], and the inverse transforms written for each variable A, 

T and A. In the inverse transform expressions, the expressions in 

a, 6 and y defining the singularities of the integrand are automatically 

the slowness surfaces. 

For the variable A, the associated slowness surface is the G = 0 
surface where G is given in eqn. (18) . This surface will be called 
SH. For the variables T and A, the remaining slowness surfaces are 
represented by H ■ 0 where H is given in eqn. (20). However, there 
are two surfaces that satisfy H = 0, resulting in two slowness 
surfaces [2,5]. These two surfaces will be called SV and P. 

The physical meaning of the slowness surfaces SH, SV and P can 
be better pursued if the isotropic medium case is considered. For 
isotropic materials 
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a 1 - a 2 , = a 5 and a 3 = a 1 - a $ (47) 

Further it can be shown that all three slowness surfaces are spherical. 
SV corresponds to purely rotational waves (purely transverse waves 
if plane waves are considered) ; SH corresponds to purely rotational 
waves (purely transverse waves if plane waves are considered) ; and P 
corresponds to purely dilatational waves (purely longitudinal waves 
if plane waves are considered). Moreover, the surfaces SH and SV 
are identical, meaning that shear waves with any polarization have 
exactly the same behavior in isotropic media. 

Formally, the identification of the surfaces SV and SH can be 
accomplished by the corresponding intersections with coordinate axes 
as follows: SH is defined as the surface that contains the shear 

mode for plane waves in the x direction or the z direction with 
particle displacements in the y direction. SV is the surface that 
contains the shear mode for plane waves in the x and z directions 
with particle displacement in the z and x directions, respectively. 
Finally, the slowness surface P is defined as the surface that 
contains the longitudinal modes in the x and z directions. The 
correspondence between wave number components (ct, 3, y) and 
geometric coordinates (x,y,z) is maintained here. 

For transversely isotropic media the rotational and dilational 
propagation modes corresponding to the slowness surfaces are not 
pure, since in general the displacement vectors are not aligned 
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with the normals of the plane front segments that constitute the 
wave front. For wave propagation in the principal directions of the 
medium though, the waves can be purely rotational or purely 
dilatational. Specifically, the plane wave propagation along a 
principal direction, the rotational mode corresponds to a shear mode 
and the dilational mode corresponds to a longitudinal mode, just 
as for the isotropic medium. The identification of the surfaces SH, 

SV and P is done as explained above, by the intersections with the 
principal directions. .An illustration of the slowness surfaces for 
a transversely isotropic medium is shown in Fig. 2. Observe that 
the values of the slownessess at the intersections of the surfaces 
with the coordinate axes are the inverses of the phase velocities 
presented in the previous section. The character of the corresponding 
plane waves along the principal directions is also sketched, showing 
the propagation and particle displacement directions. 
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GROUP VELOCITY 


The group velocity U (or velocity of energy propagation) for 
the general anisotropic medium can be written according to Rayleigh's 
energy argument [13] extended to three dimensional propagation as 
[ 1 ] 


VG 



(48) 


where G is given by eqn. (18) and U has the direction of _VG or, 
equivalently, the same direction as the normal to the slowness surface. 

For points on the slowness surface G = 0, eqn. (48) can also 
be written as 


uVG 


U = 


k-VG 


(49) 


where k is the wave number vector (a, 8, y) • Observe that the 
correspondence between the wave number vector and the slowness vector 
is [1,3] 


k = cos 


(50) 


such that k and have the same orientation. The phase velocity in 
the k direction can be represented vectorially by 



( 51 ) 
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If eqns. (49) and (51) are compared, it can be seen that 


U • k = c_ • k (52) 

Thus, the resultant of U in the £ direction is £, which is 
equivalent to saying that the velocity of energy propagation in 
the direction normal to the wave front is equal to the phase 
velocity. 

Eqns. (48) through (52) are also valid for the two other 

slowness surfaces. The same results are obtained if G is replaced 

\ 

by H in the eqns. (48) through (52), provided that Ic is taken 
accordingly. Thus, there are three phase velocities and three 
group velocities for each direction in the medium. 

Although eqn. (52) is satisfied, eqns. (49) and (51) indicate 
that U and £ are different. There is a component of energy velocity 
in the direction parallel to the wave front, meaning that energy 
transmission is not in the same direction as wave motion. This is 
true for media where the phase speed varies with direction [1]. 

Thus, plane waves oblique to the principal directions can propagate 
along their own normals only if energy is also supplied in a direction 
parallel to the wave fronts. 
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APPLICATION EXAMPLE 


Material Description 

As an application example, a unidirectional composite material 

was chosen. The material was 3M Scotchply type 1002 fiberglass epoxy 

consisting of unidirectional E-glass in a 165° C curing epoxy 

matrix. The composite prepreg tape was cured in a heated press at 

6 2 

0.69 10 N/m . The resulting resin content was 36 percent by 
weight [18] . The properties of such composites depend not only on 
the volume fraction of the components but also on the fabrication 
method. Since knowledge of the elastic properties of the material 
was necessary for the numerical calculations, samples were 
fabricated for experiments. Ultrasonic methods [14,15,16,17] were 
used for the experimental measurements, relating wave speeds and 
elastic constants c^ . The wave speeds were determined in the 
through - transmission configuration using tone bursts to generate 
either longitudinal or shear waves. The tests were performed at 
frequencies from 0.5 to 2.25 MHz in increments of 0.25 MHz. 

The test specimens were cut in the form of rectangular prisms 
with uniform square cross sections from a 25.4 mm thick composite 
plate. The cross-sectional area was square with 12.7 mm sides. 

The lengths of the specimens ranged from 2.54 mm to 25.4 mm. The 
test specimen axes were oriented at various angles with respect 
to the fiber direction. The elastic constants were determined 
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from the measured longitudinal and shear wave speeds. The 
experimental values for the elastic constants are given in Table 2. 

The convention adopted for the transversely isotropic medium 

in Fig. 1 was maintained here, so the z axis corresponded to the 

3 

fiber direction. The density was 1850 kg/m according to the 
manufacturer’s data [18]. 

Slowness Surfaces 

The slowness surfaces for the fiberglass epoxy material were 
calculated from G = 0 and H * 0, where G and H are given in eqns. 

(18) and (20) respectively. The forms of these equations indicate 
that both the surfaces are symmetric about y for a rectangular 
coordinate system a, 8, and y. Thus, in order to determine the 
surfaces, only their intersections with the a-y plane need to be 
calculated. Then, the entire slowness surfaces can be generated 
by revolving the intersection curves around the y axis. Additionally, 
the intersection curves of the slowness surfaces G = 0 and H = 0 
with the a-y plane are symmetric with respect to a axis, which can 
be checked by setting 8 = 0 in eqns. (18) and (20). By the 
combination of all symmetries, only 1/4 of the intersections of 
the surfaces on the a-y plane need to be calculated. The positive 
quadrant was chosen for the calculation. 

The ranges of the variables are limited by the values shown 
in Fig. 2 which can be obtained by setting 8=0 and alternately, 
one at a time, a = 0 or y = 0 in eqns. (18) and (20). For the positive 
quadrant the limits are: 
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For slowness surface SH: 0 <— < l/(a. 0 < — < l/(a-)^^; 

— in— 4 — a) — 5 

For slowness surface SV: 0 <— < 1/ (a )^^ 0 <— < /(a_)^^; and 

— oj — 5 — oj — 5 

For slowness surface P: 0 <— < l/(a.)^^ 0 < Y < l/(a 0 )^^. 

— to — 1 — id— 2 

The surfaces are generated by the successive combination of the 
pairs (ot,y) satisfying eqn. (18) for SH and eqn. (20) for SV and P. 

Fig. 3 shows the slowness surfaces SH, SV and P for the fiber- 
glass epoxy composite. In Fig. 3, points A and B are the so-called 
parabolic points on the slowness surface. At these points, the 
gaussian curvature based on eqn. (23) is zero. Geometrically, these 
are inflection points on the slowness surface. For the points of 
zero gaussian curvature, the amplitude coefficient X n (which corresponds 
to the inverse of a decay factor along the direction of the normal 
to the surface) goes to infinity which would give infinite 
displacement. The physical implications of the existence of such 
points can be better understood by analyzing the wave surfaces. 

From the values obtained for slownesses, the velocity surfaces 
were calculated using eqn. (44) and are shown in Fig. 4. The 
velocity surfaces, namely V(SH), V(SV) and V(P) correspond, respectively, 
to the "inverses" of the slowness surfaces SH, SV and P. 


Wave Surfaces 


The wave surfaces can be directly obtained from the slowness 

surfaces. It can be shown [1] that the coordinates x, y and z 

of the wave surfaces can be calculated from -G, /G, ; -G, /G, ; 

a o) S w 

and -G, /G, , respectively, for the wave surface corresponding to 

y (i) 
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the slowness surface SH, and by -H >a /H,^;-H,g/H,^; and -H,^/H,^ 
for the wave surfaces corresponding to the slowness surfaces SV 
and P. Fig. 5 shows the wave surfaces as calculated from the 
slowness surfaces SH, SV and P. The wave surface W (SH) , which 
corresponds to SH, has the smallest phase velocities. The wave 
surface W(P), which corresponds to P, has the largest phase 
velocities. 

As shown in Fig. 5, the wave surface corresponding to SV has 
two finite cuspidal edges [1,5]. The points at the tips of the 
cuspidal edges correspond to the inflection points A and B shown 
in Fig. 3. The directions defined by the lines from the origin 
to the tips Of the cuspidal edges (Fig. 5) represent the 
directions for which the solution is not valid. These directions 
can be found to be at 41.970° and 62.545° with respect to the z 
axis. Since the wave fronts are geometrically similar in time, 
the cuspidal edges will be always located in the same direction 
with respect to the coordinate axes . Observe that for directions 
between the angles defined above, radii from the origin to the 
W(SV) wave surface may assume three distinct values as can be seen 
in Fig. 5, indicating that there are three contributions for the 
displacements of points in this region of space. The contributions 
arise because it is possible to find three points on the 
corresponding slowness surface where the normals have the same 
direction, according to the method of stationary phase. At a 


- 30 - 



point P in a direction in space between the angles 41.970 8 and 62.545, 
three plane wave fronts can be seen traveling along the direction 
OP from the origin (where the excitation is located). Thus, it is 
possible to have plane wave fronts of different orientations and 
different phase velocities passing through the same point in the 
medium at different times. In reality these plane wave fronts 
are segments that constitute the actual wave front. 
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Displacement Calculations 


The displacements in a transversely isotropic medium subjected 
to a point force excitation as shown in Fig. 1 were also calculated 
for the fiberglass epoxy material. 

A computer program was written that calculates the displacements 
at a given point P(x,y,z) in the material. The following steps were 
performed; 

(1) Search for the point or points on each of the slowness 
surfaces in which the outward normals are parallel to the direction 
OP. This is done by comparing the direction defined by OP with the 
direction of the gradients to the slowness surfaces. The program 
retains the closest values to the direction given (one value for 
slowness surface SH; up to three values for slowness surface SV and 
up to two values for slowness surface P, in accordance with the 
possible number of points having the same normal direction) . The 
wave numbers are determined from the coordinates of the selected 
points on the slowness surfaces via eqns. (18) and (20). 

(2) Check the signs of the gaussian curvature and the 
direction of the gradient vector of the slowness surfaces, thus 
establishing the resulting phase coefficient A^. 

(3) Substitute into eqns. (36) through (38) the values found 

for a, B and v and A above and calculate the contributions of each 

n 

slowness surface to the displacements u, v and w. 
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Points in space were chosen and the corresponding displacements 
calculated. For simplicity the points were placed on the coordinate 
planes. The set of points chosen were: 

2 2 2 

(1) On the x-z plane, points along the arc x + z = 4 m , and 

2 2 2 

(2) On the x-y plane, points along the arc x + y = 4 m . 

The frequency was set at 1 MHz and a load (body force) of 
unit amplitude (F = 1 N/kg) was applied. From eqns. (15), (16) and 
(17) it can be seen that if any of the gaussian curvatures is zero, 
the corresponding displacement function assumes an infinite value. 

It is known physically that no point in space may undergo an infinite 
displacement since the amplitude of excitation F q is finite. It can 
be shown [2] mathematically that the discontinuities can be eliminated 
and the decay of the displacement amplitudes is found to vary as 
r In order to generate the amplitude coefficients, the points of 

discontinuity on the slowness surface SV are "avoided", which is 
equivalent to excluding the points A and B from the calculation set. 

Figs. 6 through 9 show the polar diagrams of the displacement 
amplitudes for the points 2 m away from the origin on the x-y plane. 
Figs. 10 through 13 show the polar diagrams of the displacements 
amplitudes for the points 2 m away from the origin on the x-z plane. 
Observe that by reasons of loading symmetry, v displacements for 
points on the x-z plane are zero as well as w displacements for points 
on x-y plane. The curves are plotted separately according to the 
slowness surface that generates the contribution. 

For points on the x-y plane, the contributions of the slowness 
surface SV are zero which implies that there are no shear waves 
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propagating along x direction with particle displacement in the z 
direction. Also, there are no shear waves propagating along the 
y direction with particle displacement in the z direction. 

The contributions from the slowness surface P, Figs. 6 and 
8, show that (1) there are pure longitudinal waves traveling along 
the x direction and this is the direction of maximum amplitude for 
u displacements (as expected since x is coincident with the force 
line) ; (2) no longitudinal waves travel in the y direction (zero 
amplitude component of v displacement for points in the y direction) , 
and no shear waves with y polarization travel in the x direction 
(v displacements are zero in x direction) . 

The contributions from the slowness surface SH, Figs. 7 and 9, 
show that (1) there are pure shear waves in the y direction (with 
x polarization direction) but since v displacements along the y 
direction are zero, these are standing waves (y direction vibrates 
as a string) ; (2) no longitudinal waves travel in the y direction 
(zero amplitude component of v displacements for points in the 
y direction) , and no shear waves with y polarization travel in the x 
direction (v displacements are zero in x direction) . 

For points on the x-z plane, the contributions of slowness 
surface SH are zero which implies that there are no shear waves 
propagating along the x direction with particle displacement in 
the y direction. Also there are no waves propagating along the z 
direction with particle motion in the y direction. 
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The contributions from the slowness surface SV, Figs. 10 and 
11, show that (1) there are pure shear waves in the z direction (with 
x polarization direction) but since the w displacements along the 
z direction are zero these are standing waves (z direction vibrates as 
a string) ; (2) there are no shear waves traveling along the 
x direction (zero amplitude component of w displacement in the x 
direction) . 

The contributions from the slowness surface P, Figs. 11 and 13, 
show that (1) there are pure longitudinal waves propagating along 
the x direction and the u displacement amplitude along the x direction 
is larger than the corresponding u displacement amplitude along the 
z direction; (2) there are pure shear waves in the z direction (with 
x polarization direction) but since the w displacements along the 
z direction are zero these are standing waves. 

By observing the polar diagrams of displacements (Figs. (6) 
through (13)), it can be seen that there are characteristic directions 
defining maxima for the displacement amplitudes. For comparison 
purposes these directions are listed in Table 3. It can be observed 
that there are two maxima for the contributions of slowness surface 
SV corresponding to the inflection points A and B in Fig. 3. For 
the contributions of the slowness surface P, the maxima occur at 
two other distinct positions. For the points on the plane x-y, 
the maxima occur along the line of force (x axis itself) for u 
displacements except for the small contribution of SH (see Fig. 7), 
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and at 45° for v displacements as expected, since x-y is isotropic. 
Slowness furfaces not listed in Table 3 have zero contributions to 
the corresponding displacements. 

Recalling the relationship between strain energy and displacement 
amplitude (energy is proportional to the square of amplitude) , it 
is clear that for points in regions of larger amplitudes the energy 
will assume the higher values. The difference between anisotropic 
and isotropic media is that the energy density for anisotropic media 
is not a uniform function for equal solid angles having vertexes at 
the origin, constructed around different directions in space. 

Lighthill [2] referred to this phenomenon saying that the energy is 
"confined" to some cone in space, or in other words, most of the 
energy travels along preferential directions. The preferential 
directions for the unidirectional fiberglass epoxy composite are 
shown in Table 3. It is understood that the maxima are the center 
directions of the regions around which the energy is confined. When 
the numerical process is carried out the preferential directions are 
automatically determined. 
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CONCLUSIONS 


The far-field displacement pattern in an infinite transversely 
isotropic medium subjected to an oscillatory point force was determined 
and evaluated for a specific fiberglass epoxy composite. The solution 
describes the stress wave field in terms of the geometric aspects 
of the disturbance spreading by means of slowness and wave surfaces. 

The solution for displacements allows the prediction of the amplitude 
distribution as shown by the construction of polar diagrams. 

It was seen that the energy from a point source travels along 
radial directions even if the wave surfaces assume complicated shapes, 
meaning that the phase and group velocities are different. 

The existence of preferential directions is an important aspect 
to be considered when experimental tests are to be designed. 

Knowledge of the displacement field allows a better choice to be 
made for the positioning of components of the measuring system. 

There appear to be no restrictions regarding the applicability 
of this method to other types of anisotropy. If appropriate symmetry 
relations can be applied as in the case of eqn. (7) , the same 
procedurescan be followed for the displacement field solution. This 
possibility is attractive for use in the description of the stress 
wave field of filamentary composite materials with fiber arrangements 
other than unidirectional. 

Another possible application is the study of acoustic emission 
(AE) phenomena in composites or other anisotropic materials. It is 
known that acoustic emission is a transient disturbance of relatively 
small dimensions occurring inside the material, and normally 
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containing a range of frequency components. If the method is applied 
for each component of frequency, superposition can be used for the 
determination of the field resulting from the AE source. Knowledge 
of the displacement field might allow inferences regarding the 
source such as strength and orientation, and consequently perhaps 
the degree of damage within the material. 
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TABLE 1 Velocities for Plane Waves along Principal 
Directions of Transversely Isotropic Media, 
(x-y is isotropic plane) 


Direction of 
Propagation 

Direction of Displacement Vector 
x y z 

X 

(a/' 2 

u 4 > 1/2 

(a 5 ) 1/2 

y 

(V 1/2 

(a/' 2 

(a 5 ) 1 « 

z 

u 5 ) 1/2 

(a/ /2 

(a/' 2 


The a^ are defined in eqn. (13). 
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TABLE 2 Elastic Constants for Unidirectional Fiberglass 
Epoxy Composite (Scotchply 1002) . 


Elastic 

Constant 

Computed Value from 
Wave Speed Measurements 

C 11 

10.580 

10 9 

N/m 2 

C 12 

4.098 

10 9 

2 

N/m 

C 13 

4.679 

10 9 

N/m 2 

C 33 

40.741 

10 9 

N/m 2 

C 44 

4.422 

10 9 

N/m 2 
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TABLE 3 Points of Maximum Displacement Amplitude for Lines 
on x-z Plane and x-y Plane 


Line Equation 

Displm 

Slowness 

Surfaces 

Maximum 

Location 

Order of 
Magnitude 

of Amplitudes (***) 

2 2 . 
x + z =4 

u 

SV 

41°; 62.5° (*) 

io -8 

P 

90° (*) 

10 -8 

w 

SV 

42°; 62.6° (*) 

io -9 

P 

15 0 (*) 

io " 9 

2 3 . 

x + y = 4 

u 

SH 

0 ° (**) 

io " 8 

P 

90° (**) 

10 -9 

V 

SH 

45° (**) 

io -9 

P 

45 ° (**) 

10“ 9 


(*) Angles with respect to z axis (in degrees) 

(**) Angles with respect to y axis (in degrees) 

(***) These amplitudes correspond to a unit body force. 
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a- wave number vector component in x direction 
y- wave number vector component in z direction 

— *■ Propagation direction of associated plane wave (along 
principal directions) 

— ♦ Particle displacement of associated plane wave (along 
principal directions) 


Fig. 2 Plane representation of slowness surfaces for transversely 
isotropic medium where z axis corresponds to the larger 
elastic constant value and x-y is the isotropic plane. 
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Component of slowness in z direction y/cu(s/m) 



Component of slowness in x direction a/cu (s/m) 

Fi^. 3 Slowness surfaces SH, SV and P for unidirectional 
fiberglass composite for positive x-z quadrant. 
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Component of phase velocity in z direction c 2 (m/s) 



Component of phase velocity in x direction c x (m/s) 


Fig. 4 Velocity surfaces V(SH), V(SV) and V(P) for 
unidirectional fiberglass composite for 
positive x-z quadrant. 
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Distance traveled by wavefront in z direction 
in I second z(m) 



Distance traveled by wavefront in x direction 

in I second x(m) 


Fig. 5 Wave surfaces W(SH), W(SV) and W(P) for unidirectional 
fiberglass composite for positive x-z quadrant. 
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Amplitude component in y direction (m) 



Amplitude component in y direction (m) 









Amplitude component in y direction (m) 



Fig. 8 Polar diagram for v displacement amplitudes for 
points along the line x 2 + y 2 = 4m 2 , due to 
slowness surface P only. 


- 51 - 


Amplitude component in y direction (m) 



Amplitude component in z direction (m) 


2.8 x 10“ 


0 0.4 0.8 1.2 1.6 2.0 2.4 2.8x10 

Amplitude component in x direction (m) 


Fig. 10 Polar diagram for u displacement amplitudes for 

points along the line x 2 + z 2 = 4m 2 , due to slowness 
surface SV only. 


Amplitude component in z direction (m) 




Amplitude component in z direction (m) 


6XICT 9 



Amplitude component in x direction (m) 


Fig. 12 Polar diagram for w displacement amplitudes for 

points along the line x 2 + z 2 = 4m 2 , due to slowness 
surface SV only. 
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Amplitude component in z direction (m) 
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